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SARLAC:  A  RELATIVISTIC  ELECTRON  BEAM  CODE 

We  have  recently  written  a  number  cf  simulation  codes  to  test  various 

k 

aspects  of  the  resistive  hose  instability  in  high  energy  electron  beams 
propagating  in  resistive  plasmas.  Most  methods  used  previously  for 
treating  the  instability  were  restricted  to  small  instability  amplitudes. 
These  are  considered  to  be  of  practical  interest  because  large  amplitude 
hose-like  oscillations  quickly  destroy  the  integrity  of  the  the  beam,  and 
because,  under  appropriate  circumstances,  the  instability  "saturates"  in 
the  linear  regime.  That  is  to  say,  the  hose  instability  is  convective  in 
the  beam  frame,  and  therefore,  at  any  given  point  in  the  beam  may  reach  a 
maximum  which  is  still  small  followed  by  a  decay  of  the  instability  as  the 
disturbance  convects  past.  Near  term  experiments,  however,  are  frequently 
more  unstable  than  can  be  treated  by  linear  models,  so  we  have  developed  a 

particle  simulation  code  which  can  follow  the  evolution  of  an  instability 

1-3  4 

into  the  nonlinear  regime  .  Similar  codes  have  been  written  by  Godfrey 

and  Freeman"*. 

The  nonlinear  code  has  borrowed  heavily  from  two  of  our  previous 
codes,  SIMMO^  and  SIMMl^  which  are  particle  simulation  codes  for  axi- 
symmetric  beams  and  beams  with  small  amplitude  hose  motions.  The  particle 
dynamics  of  those  two  codes  are  followed  in  Cartesian  coordinates  so  the 
SARLAC  code  differs  from  these  primarily  in  the  calculation  of  the 
electromagnetic  fields.  We  have  developed  a  fast  iterative  field  solver 
which  allows  us  to  include  a  large  number  of  Fourier  modes  in  the  azimuthal 
direction. 

The  code  employs  many  of  the  approximations  found  in  most  linearized 
8-10 

propagation  models  .  The  variables  z  and  t  are  replaced  by  z  and 
C  -  ct  -  z  (the  distance  from  the  beam  head),  and  all  particles  remain  at 
constant  f,  since  v^  is  assumed  to  be  the  velocity  of  light.  The  frozen 
approximation  is  used  in  the  field  equations,  and  the  same  conductivity 
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model  used  in  the  VIPER  code  is  employed.  Beam  dynamics  are  treated  using 
standard  particle  simulation  techniques.  Current  densities,  fields,  and 
conductivity  are  calculated  on  a  polar  grid  (u,0, C)  with  u  =  /r  as  the 
radial  variable.  The  lay-down  scheme  for  the  particles  is  quadratic  in  the 
radial  and  azimuthal  variables  and  nearest  grid  point  in  the  axial 
variable. 

The  ultra-relativistic  approximations  used  in  SARLAC  lead  to  a  code 

structure  which  is  substantially  different  from  "conventional"  particle 

simulations.  Information  can  only  flow  in  one  direction;  toward  larger  C» 

Also,  since  individual  particles  always  remain  at  the  same  axial  position 

within  the  beam,  the  simulation  can  treat  one  slice  at  a  time,  thus 

reducing  the  number  of  particles  in  the  simulation  at  any  one  time  to  ~10^. 

Each  beam  slice  is  propagated  forward  in  z  until  the  maximum  propagation 

range,  z  is  reached.  At  this  point,  particles  are  loaded  into  the  next 
max 

slice,  and  the  process  is  repeated.  The  current  density  J,  conductivity  a, 
and  potentials  A  and  4>  from  the  previous  slice  must  be  read  from  disk.  The 
axial  step  AC  is  variable,  and  the  code  has  the  option  of  subgridding  the 
field  and  conductivity  integrations  on  a  finer  axial  mesh  than  is  used  for 
the  particles.  All  diagnostics  are  done  with  post-processors.  The 

7  8 

dimensionless  units  used  in  VIPER  and  SIMM1  are  employed  throughout  ’  . 


ELECTROMAGNETIC  FIELDS 


The  frozen  approximation  to  Maxwell's  equations  is  performed  in  a 
11 

gauge  suggested  by  Lee  .  The  equations  are 
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with  the  frozen  condition 


d  -  _  da  _  d(f>  _  d  A  ± 


dz 


dz 


dz 


(3) 


and  <t  =  A  -  <}>. 
z 


The  boundary  conditions  for  a  and  <|>  are 
a(r)  =  0  ,  C  =  0 

<j>(r)  =  0  ,  C  =  0 

®<0  =  0  ,  r  =  R 


max 


<KO  =  0  ,  r  =  R 


max 


These  conditions  correspond  to  a  beam  propagating  at  the  speed  of  light  in 
a  perfect  conductor  of  radius  R 


max 
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The  equations  are  similar  in  form  to  the  EMPULSE  equations  with  an 
2  2 

additional  term  9  #/9£  in  the  first  equation.  A  fully  implicit  method  for 

13 

solving  these  equations  has  been  developed  by  Hui  .  That  field  solver 
Fourier  analyzes  the  azimuthal  dependence  of  all  quantities  into  a  series 
of  modes  exp(im0)  and  performs  a  full  complex  matrix  inversion,  which  is 
extremely  time  consuming  and  thus  impractical  for  long  simulation  studies. 
The  major  advance  of  the  SARLAC  code  is  the  development  of  a  field  solver 
which  does  not  require  a  complete  matrix  inversion.  The  SARLAC  field 
solver  uses  a  predictor-corrector  method  which  iterates  about  a  solution 
obtained  by  assuming  that  the  axi-symmetric  (m  =  0)  conductivity  dominates 

■  t 

the  solution.  The  m  =  0  mode  of  any  positive  definite  function  is  always 
larger  than  any  other  single  mode  and  in  the  case  of  beam  generated 
conductivity  which  is  generated  all  along  the  beam  axis,  this  mode  is  large 
compared  to  the  other  modes  even  for  large  excursions  of  the  beam  from  axi- 
symmetry,  as  long  as  the  front  of  the  beam  is  on  the  axis. 

Consider  the  first  equation,  and  write  it  in  the  form 


if-*-- 


(4) 
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Integrate  this  equation  over  the  interval  C  <  C  <  C  , 

n  --  ^n+1 


to  obtain 


n+1  In  n-1  1  (.  -ctAC 

«  =  \a  ~  &  oTTF  1  -  e  S 


') ♦  Jt  h  (»« -  I1  -  e'"AC)) 


i(a  +  ♦>  7  (•«  - 11  -  •"*))}/(»  -  fa  t1  -  -“aac)).  (5> 


where  the  superscript  n  represents  the  function  evaluated  at  C  .  The 
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conductivity,  a,  is  evaluated  in  the  interval  (^n>^n+^)*  Rewrite  Eq.  (5) 
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Note  that 


AC/ a  for  ct  large 


— >  ( AC)  for  a  small , 


Equation  (2)  is  simply  differenced  to  give 


n  n+1  n 
,2  a  -a 


-  \  — AC^-=  -  V<°V> 


Note  that  we  have  omitted  the  superscripts  from  some  of  the  a  and  <{>  terms 
in  Eqs.  (6)  and  (9).  By  choosing  these  terms  at  the  nth  or  (n+l)th  (or 
some  combination  of  these  levels),  the  differencing  can  be  made  explicit, 
or  implicit  to  some  formal  accuracy.  Our  experience  has  been  that  the 
algorithms  for  these  equations  are  numerically  unstable  if  they  are 
explicitly  differenced.  An  implicit  differencing  can  eliminate  the 
instability,  but  at  the  expense  of  a  complicated  matrix  inversion  due  to 
the  azimuthal  coupling  of  a  with  the  potentials.  To  avoid  this,  we  have 
chosen  instead  to  to  rewrite  Eqs.  (6)  and  (9)  as 

an+1  -  /„(*)  =  F(a\an-'-,Jb,4>) 

+  v](«  +  40  [/(</>)  -  /oW)]  (10) 

and 


-  <5  — -  V (•<>*+) -  -  V<*  -  W’  (11) 

where 


fo  to  =  “  j0*f  (°(6)d  6,  a0  =  ■£-  f  a  (0)  dd 


(12) 


Equations  (10)  and  (11)  are  formally  the  same  as  Eqs.  (6)  and  (9). 
However,  since  cro  does  not  vary  in  azimuth,  the  left-hand  sides  of  these 
equations  can  be  evaluated  at  the  upper  level  without  involving  convolution 
sums,  which  leads  to  a  tri-diagonal  form  for  #.n+^  and  4>n+^" .  The  right-hand 
sides  can  be  evaluated  explicitly  since  the  functions  are  known.  The 
simplest  differencing  scheme  leads  to  a  first  order  algorithm.  We  have 
chosen  a  predictor-corrector  method  which  is  accurate  to  second  order  and 
quite  stable.  For  the  sake  of  brevity,  we  will  not  go  into  the  details  of 
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the  differencing,  but  they  are  easily  reproduced.  We  solve  the  equations 

in  the  Fourier  transformed  space  (u,m,C)  always  keeping  the  right-hand 

sides  explicit.  The  right-hand  sides,  then,  can  be  evaluated  in  (u,0, C) 

thus  avoiding  convolution  sums  on  the  right-hand  sides  as  well.  The  method 

of  adding  and  subtracting  averaged  terms  to  gain  stability  without  full 

14 

matrix  inversions  is  similar  to  that  used  by  Harned  for  a  different  set 
of  equations. 

NUMERICAL  ISSUES 

In  SARLAC,  the  number  of  modes  Nq  and  radial  mesh  size  Au  remain  fixed 

1/2 

throughout  a  run.  Typically,  Nq  =  16  or  32,  and  Au  =  0.125  aQ  ,  where 

a  is  a  characteristic  initial  beam  radius.  The  axial  grid  spacing  AC  is 
o 

specified  for  each  slice  at  the  beginning  of  the  run.  In  general,  AC  is 

allowed  to  increase  with  C  since  (at  least  in  the  linear  regime)  the  C- 

2  2 

variation  is  characterized  by  the  dipole  decay  length,  na(r=0)aQ  /2c  , 
which  usually  increases  monotonically  throughout  the  pulse.  However,  field 
solver  tests  have  shown  that  the  axial  step  size  must  often  be  reduced  when 
the  beam  displacement  is  large.  This  is  accomplished  by  subgridding  the 
field  and  conductivity  integrations.  In  most  cases,  AC  is  chosen  to  be 
small  enough  that  subgridding  is  rarely  involved. 

The  beam  current  density  is  intrinsically  noisy  because  of  the 
statistical  fluctuations  arising  from  the  small  number  of  particles  in  each 
u-0  grid  cell.  This  is  particularly  troublesome  near  u  =  0.  Increasing 
the  number  of  simulation  particles  per  slice  reduces  noise  problems  but  is 
computationally  expensive.  Other  methods  which  we  have  employed  include 
accumulating  current  densities  on  a  coarser  radial  mesh  than  is  used  for 
the  field  solver  and  interpolating,  averaging  over  the  first  few  radial 
grid  points,  and  using  an  azimuthal  filtering  technique  near  the  origin. 
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Originally,  we  assigned  random  initial  values  of  0  to  the  particles  but 
found  that  this  procedure  resulted  in  large  initial  noise  levels  for  the 
hose  instability  and  in  substantial  drifts  in  the  beam  head.  The  noise 
effects  are  reduced  by  loading  the  particles  in  pairs  on  opposite  sides  of 
the  beam.  If  the  velocities  are  also  loaded  symmetrically,  the  m  =  1 
azimuthal  mode  is  eliminated  in  the  initial  stages  of  propagation.  A  small 
specified  perturbation  can  be  added  to  all  particles  in  a  given  beam  slice 
to  start  the  hose  instability  in  a  controlled  manner.  Higher  order  Fourier 
modes  can  be  suppressed  by  loading  four  or  more  particles  with  the 
appropriate  symmetry.  The  elimination  of  higher  order  modes  has  not  proven 
particularly  useful  since  the  nonlinear  coupling  of  these  modes  is  usually 
too  weak  to  introduce  significant  hose  growth. 

The  scattering  of  beam  electrons  by  the  neutral  gas  is  known  to  play  an 
important  role  in  the  evolution  of  the  beam.  SARLAC  uses  an  algorithm 
originally  developed  by  Chambers6’15  and  modified  by  Hughes  and  Godfrey16 
to  provide  a  more  accurate  representation  of  the  scattering  process.  Each 
beam  particle  is  periodically  scattered  through  a  randomly  chosen  angle 
whose  characteristic  magnitude  is  determined  by  the  energy  and  the  gas  density. 
After  an  initial  transient  phase,  the  beam  reaches  a  quasi-static 
equilibrium.  The  beam  radius  then  expands  slowly  due  to  scattering.  If 
beam  particles  are  loaded  in  pairs,  a  straightforward  application  of  the 
scattering  algorithm  will  eventually  introduce  significant  noise  and  drifts 
at  the  beam  head.  These  effects  can  be  eliminated  by  scattering  the 
particles  in  pairs.  The  random  velocity  Av.  applied  to  a  given  particle  at 
a  given  z  step  is  balanced  by  adding  -Av^  to  the  particle  with  which  it  was 
originally  paired.  This  technique  has  been  highly  successful  in  practice. 
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The  simulation  code  is  best  suited  for  treating  cases  in  which  the  hose 
displacement  is  a  few  times  the  nominal  radius  aQ.  The  coordinate  system 
is  chosen  to  have  the  finest  resolution  in  both  the  radial  and  azimuthal 
directions  near  the  origin.  For  extremely  large  beam  oscillations,  when 
the  beam  displacement  reaches  a  substantial  fraction  of  the  wall  radius, 
the  accuracy  of  the  simulation  is  reduced,  and  the  field  solver  is 
sometimes  subject  to  numerical  instabilities.  The  field  solver 
instabilities  appear  to  be  triggered  by  conditions  in  which  the  local 
conductivity  centroid  gets  far  enough  off  the  coordinate  system  axis  that 
the  conductivity  is  not  dominated  by  the  m  =  0  mode.  Usually  this 
conductivity  is  generated  by  avalanche  due  to  strong,  localized,  electric 
fields;  such  fields  can  arise  when  the  hose  motion  is  quite  nonlinear1. 
Evidence  for  very  strong  electric  fields  associated  with  nonlinear  hose 
motion  has  been  seen  in  the  ETA  experiments17,  so  the  strong  fields  may  be 
physical  (up  to  a  point).  Considerable  effort  has  been  made  to  make  the 
field  solver  more  robust,  and  with  careful  differencing  we  have  had  some 
success.  We  have  also  found  that  these  problems  can  be  mitigated  by  using 
small  C.  grids  in  regions  where  there  are  large  hose  amplitudes.  Even  so, 
we  believe  the  code  to  be  best  suited  for  moderate  hose  oscillations. 


NUMERICAL  RESULTS 


We  have  run  the  code  under  a  variety  of  conditions.  We  show  here  the 

results  of  two  runs;  one  for  small  perturbations  in  which  the  hose  stays  in 

the  linear  regime,  saturates  and  decays,  and  one  with  moderate  initial 

perturbations  for  which  the  hose  grows  and  becomes  nonlinear.  The 

parameters  for  both  sets  of  runs  are 

a  =  .5  cm,  the  beam  radius, 
o 

I  =10  kA,  the  beam  current. 

Y  =  100,  the  beam  energy. 

=  15  cm,  the  beam  current  rise  length. 

a  =  81  a  ,  the  outer  radius  of  the  simulation, 
w  o 

y  =  y  sin  211  ((£-£  )/C  )*  the  initial  perturbation  over  the 
•'pert  o  o  o 

range 

C  <  C  <  1.5  C  ,  C  =  10  cm. 

0  o  o 

(Note,  this  perturbation  is  in  the  y  direction.) 

For  the  first  run  we  used  a  very  small  initial  hose  amplitude  y  =10  \ 

oo 

which  kept  the  hose  oscillations  linear  over  the  length  of  the  beam.  For 

this  case,  the  hose  instability  grows  and  saturates  as  seen  in  Fig.  1. 

Figure  2  shows  a  comparison  of  the  saturation  amplitude  at  various 

distances  from  the  beam  head  with  the  results  of  the  linearized  VIPER  code, 

18 

which  uses  the  multi-component  model  to  represent  the  particle  dynamics 

approximately.  The  oscillation  frequency,  growth  rate,  and  saturation 

8 

amplitude  agree  quite  well  with  the  VIPER  code  . 

The  parameters  chosen  for  the  second  run  were  the  same  except  for  a 

_2 

much  larger  initial  amplitude  yQ  =  10  ao,  so  that  the  hose  oscillations 
would  become  nonlinear.  Figure  3  shows  the  growth  of  the  hose  through  the 
x  and  y  centroids  of  the  beam.  The  dashed  line  is  the  y  centroid  which  is 


much  larger  than  x  because  the  perturbation  is  initialized  in  y.  After  the 


hose  displacements  reach  the  order  of  the  beam  radius,  the  frequency  of  the 


-  oscillation  decreases.  This  is  because  the  beam  is  spreading  in  radius  and 

. 

i 

the  wavelength  of  the  oscillation  scales  as  the  radius.  Figures  4-6  are 


beam  particle  plots  at  various  distances  of  propagation.  We  can  see  the 


development  of  the  hose  instability  and  the  loss  of  the  beam  pinch  as  the 


instability  becomes  large. 

I 

I 
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Figure  2.  Plot  of  the  saturated  hose  amplitude  as  a  function  of  the 

distance  back  from  the  point  of  the  initial  perturbation.  The 
SARLAC  results  are  marked  with  *'s,  while  the  VIPER  results  are 


marked  with  +  's. 
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Figure  3.  Plot  of  the  beam  x  and  y  centroids  at  C  =  50  cm  as  a  function  of 

the  beam  propagation  distance.  The  initial  perturbation  is 
-2 

yQ  =  10  so  that  the  hose  instability  becomes  nonlinear.  When 
the  hose  motion  becomes  large,  the  beam  radius  increases  and  the 
frequency  of  the  oscillation  decreases.  Again,  the  dashed  line 


ng  C  but  this  is  not  indicated  by  increasing  numbers 


e  particles  are  weighted  with  the  charge. 


cm 


X-pos  vs  Zeta 


Z  =  600.  Zeta  (cm) 

Figure  6.  Plot  of  the  beam  particle  positions  at  z  =  600  cm.  The  hose 


perturbation  has  grown  so  large  that  the  tail  of  the  beam  is 
completely  disrupted. 
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